Breareal

Analyser av breareal

Importerer data

Det er fire datasett - breatlaset, 50, regioninndelingen, og omriss av Norge. Vi trenger ikke fjellmasker da vi kan anta at alle isebreer ligger i fjellet.

Breatlas

Her er det nye breatlasset med breareal fra 2018 og 2019 hentet fra Sentinell (ref Liss Marie Andreassen, NVE)

breatlas <- sf::st_read('/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/breatlas/breatlas_2018_2019/Breatlas_20182019_temp20210922_lma.shp')

## Reading layer `Breatlas_20182019_temp20210922_lma' from data source 
##   `/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/breatlas/breatlas_2018_2019/Breatlas_20182019_temp20210922_lma.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6915 features and 27 fields (with 4 geometries empty)
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -9874.591 ymin: 6624398 xmax: 810132 ymax: 7836994
## Projected CRS: ETRS89 / UTM zone 33N

CRS: UTM 33N ETRS89 Det er 6915 polygoner

plot(breatlas[1], main = "Breatlas 2018-2019", col = "black")

Arealet ser større ut en det er fordi det er lagt på kantlinjer.

breatlas$area <- st_area(breatlas)
par(mfrow=c(1,2))
hist(breatlas$area)
plot(breatlas$area)

De fleste polygonene er små. Og så er det noen få veldig store. Den største er 50 km2.

Kart over Norge

nor <- readRDS('../data/norway_outline.RDS')%>%
  st_as_sf()%>%
  st_transform(crs=crs(breatlas))
plot(nor$geometry, axes=T, main = "Breatlas 2018-2019")
  plot(breatlas$geometry, add=T, border = "blue")

UTM33 WGS84

N50

Her er utstrekningen til breer i perioden 1952 til 1985, basert på gamle kart.

n50 <- sf::st_read('/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/breatlas/n50/cryoclim_GAO_NO_1952_1985_UTM_33N.shp')%>%
  st_transform(crs = crs(breatlas))

## Reading layer `cryoclim_GAO_NO_1952_1985_UTM_33N' from data source 
##   `/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/breatlas/n50/cryoclim_GAO_NO_1952_1985_UTM_33N.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 7141 features and 30 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -7278.562 ymin: 6623024 xmax: 812207.8 ymax: 7841654
## Projected CRS: WGS 84 / UTM zone 33N

plot(nor$geometry, axes=T, main = "n50 - 1952-1985")
  plot(n50$geometry, add=T, border = "red")

La oss plotte alt sammen oppå hverandre

myExt <- raster::extent(c(0, 100000, 6840000, 6900000))
myExt2 <- raster::extent(c(5000, 10000, 6880000, 6885000))

par(mfrow=c(1,3))
plot(nor$geometry, axes=T)
    plot(n50$geometry,  border = "orange", col = scales::alpha("orange", 1), add=T)
    plot(myExt, add=T, lwd=2)

plot(nor$geometry, xlim=c(0, 100000),
          ylim=c(6840000, 6900000),
          axes=T)
    plot(n50$geometry, add=T, col = "grey")
    plot(myExt2, add=T, lwd=3, col="orange")
    
plot(nor$geometry, xlim=c(5000, 10000),
          ylim=c(6880000, 6885000),
     axes=T)
    plot(n50$geometry, add=T, col = "red")
    plot(breatlas$geometry, add=T, col = "grey")

De røde områdene i kartet til høyre er hvor isen har trekt seg tilbake. Nå må vi dele kartlage inn etter regioner for så å sammenligne arealet i n50 med breatlaset. Hvordan vi gjør det med kalkulering av usikkerhet vet jeg fortsatt ikke.

Regioner

Henter shp med regionene

reg <- st_read("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/FINAL/Raw_data/Geografisk_oppdeling/regioner_2010/regNorway_wgs84 - MERGED.shp")%>%
  st_transform(crs = crs(breatlas))

## Reading layer `regNorway_wgs84 - MERGED' from data source 
##   `/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/FINAL/Raw_data/Geografisk_oppdeling/regioner_2010/regNorway_wgs84 - MERGED.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 4.087524 ymin: 57.75901 xmax: 31.76159 ymax: 71.38488
## Geodetic CRS:  WGS 84

plot(nor$geometry, axes=T)
  plot(reg$geometry, add=T, border = "black", 
       col = scales::alpha(c("blue", 
                             "red", 
                             "green",
                             "yellow",
                             "brown"), .2))

Før vi henter brearealet inn i tabellen for regionene må jeg prøve å fikse et problem med noen brepolygoner som krysser seg selv

table(st_is_valid(breatlas))

## 
## FALSE  TRUE 
##    68  6847

st_is_valid(breatlas, reason=T)[1:10]

##  [1] "Valid Geometry"                                  
##  [2] "Valid Geometry"                                  
##  [3] "Valid Geometry"                                  
##  [4] "Valid Geometry"                                  
##  [5] "Valid Geometry"                                  
##  [6] "Valid Geometry"                                  
##  [7] "Valid Geometry"                                  
##  [8] "Valid Geometry"                                  
##  [9] "Ring Self-intersection[100236.4931 6912507.5864]"
## [10] "Valid Geometry"

# krever lwgeom
breatlas2 <- st_make_valid(breatlas)
table(st_is_valid(breatlas2))

## 
## TRUE 
## 6915

# samme resultat:
#breatlas2 <- st_buffer(breatlas, 0.0)
#table(st_is_valid(breatlas2))

plot(nor$geometry, xlim=c(5000, 10000),
          ylim=c(6880000, 6885000),
     axes=T)
    plot(breatlas2$geometry, add=T, col = scales::alpha("blue",0.5))
    plot(breatlas$geometry, add=T, col = scales::alpha("yellow",0.5))

Dette ser ut til å ha funket fint. Denne litt kunstige inndelingen av tilgrensende polygoner er forresten bare i det nye breatlaset, ikke i n50.

breatlas <- breatlas2
rm(breatlas2)

unique(reg$region)

## [1] "Nord-Norge"      "Midt-Norge"      "Ã\u0098stlandet" "Vestlandet"     
## [5] "Sørlandet"

Fikser øæå

reg$region[reg$region=="Ã\u0098stlandet"] <- "Østlandet"
reg$region[reg$region=="Sørlandet"] <- "Sørlandet"

Breareal per region

Dagens breer

Her er en metode for å finne breareal per region.

#brealtas_reg <- st_intersection(breatlas, reg)
#saveRDS(brealtas_reg, "../data/brealtas_reg_helperfile.rds")
brealtas_reg <- readRDS("../data/brealtas_reg_helperfile.rds")

Skjekker hva som skjer med breer som ligger midt mellom to regioner

plot(brealtas_reg$geometry[brealtas_reg$region=="Midt-Norge"], 
    col = scales::alpha("red",0.5), border=NA, axes=T, 
     ylim=c(6900000, 6905000),
     xlim=c(85000, 100000))
plot(brealtas_reg$geometry[brealtas_reg$region=="Vestlandet"], 
     col = scales::alpha("green",0.5), border=NA, add=T)
  plot(nor$geometry, add=T)

Det deler polygonene ved grensen. Det er bra

brealtas_reg$area_crop <- st_area(brealtas_reg)

bretabell <- tapply(
       brealtas_reg$area_crop,
       brealtas_reg$region,
       FUN = sum)
bretabell <- bretabell/1000000 
barplot(
  bretabell,
  ylab="Breareal (km2)"
)

### Samme for N50

n50x <- st_make_valid(n50)
#n50_reg <- st_intersection(n50x, reg)
#saveRDS(n50_reg, "../data/n50_helperfile.rds")
n50_reg <- readRDS("../data/n50_helperfile.rds")

plot(n50_reg$geometry[n50_reg$region=="Midt-Norge"], 
    col = scales::alpha("red",0.5), border=NA, axes=T, 
     ylim=c(6900000, 6905000),
     xlim=c(85000, 100000))
plot(n50_reg$geometry[n50_reg$region=="Vestlandet"], 
     col = scales::alpha("green",0.5), border=NA, add=T)
  plot(nor$geometry, add=T)

n50_reg$area_crop <- st_area(n50_reg)

n50tabell <- tapply(
       n50_reg$area_crop,
       n50_reg$region,
       FUN = sum)
n50tabell <- n50tabell/1000000 


barplot(
  n50tabell,
  ylab="Breareal (km2)",
  col="grey"
)
barplot(
  bretabell,
  ylab="Breareal (km2)",
  add=T,
  col="grey30"
)

I figuren over er brearealet i referanseperioden i lyegrått og brearealet i 2018-2019 i mørkegrått. Høyden på det lysegråe partiet viser med andre ord reduksjonen i breareal i denne perioden.

Skalerte verdier

myTbl <- tibble(region          = names(bretabell),
                indikator        = bretabell,
                referanseverdi   = n50tabell)
myTbl$skalert_indikator <- myTbl$indikator/myTbl$referanseverdi

Skalerer og klipper regionene etter norgeskartet

reg$skalert_indikator <- myTbl$skalert_indikator[match(reg$region, myTbl$region)]
reg_clipped <- st_intersection(reg, nor)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

reg_clipped$skalert_indikator_r <- round(reg_clipped$skalert_indikator, 2)
tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(reg_clipped) + 
  tm_polygons(col="skalert_indikator", border.col = "white")+
  tm_text("skalert_indikator_r", size=1.5)+
  tm_shape(nor)+
  tm_polygons(alpha = 0,border.col = "black")

Det er 4 måter å kalkulere den nasjonale indikatorverdien på:

1. ta gjennomsnitt av de skalerte regionale indikatorverdiene og la alle regoinene telle likt (regional skaleringsmetode)

(alt1 <- mean(reg_clipped$skalert_indikator))

## [1] 0.5445936

2. summer indikatorverdien for hele landet og del på den summerte referanseverdien for hele landet (nasjonal metode)

(alt2 <- sum(bretabell)/sum(n50tabell))

## [1] 0.7003536

3. summere skalerte regionale indikatorverdier etter å ha gitt de ulike vekt grunnet totalarealet til regionen (arealveid metode)

reg_clipped$areal <- st_area(reg_clipped$geometry)/1000000

totalArea <- sum(reg_clipped$areal)

reg_clipped$skalert_indikator_v <- reg_clipped$skalert_indikator*
  reg_clipped$areal


sumSkalerteIndikatorer <- sum(reg_clipped$skalert_indikator_v)

(alt3 <- sumSkalerteIndikatorer/totalArea)

## 0.5884451 [1]

4. summere skalerte regionale indikatorverdier etter å ha gitt de ulike vekt grunnet fjellareal (fjellarealveid metode)

For dette trenger jeg fjellareal per region:

(fjellareal <- readRDS("../data/fjellareal.rds"))

##       Region Fjellareal
## 1 Nord-Norge   59822.74
## 2 Midt-Norge   18294.64
## 3  Østlandet   19004.18
## 4 Vestlandet   20329.63
## 5  Sørlandet    7085.62

reg_clipped$fjellareal <- fjellareal$Fjellareal[match(reg_clipped$region, fjellareal$Region)]
reg_clipped$skalert_indikator_v2 <- reg_clipped$skalert_indikator*
  reg_clipped$fjellareal
(alt4 <- sum(reg_clipped$skalert_indikator_v2)/sum(reg_clipped$fjellareal))

## [1] 0.619549

Her er alternativene oppsummert:

1. Indeksverdi 0.54. Aggregeringsmetoden er tvilsom.

2. Indeksverdi 0.7. Aggregeringsmetoden sier oss hvordan den TOTALE endringen i breareal har vært for Norge sett under ett.

3. Indeksverdi 0.59. Aggregeringsmetoden er igjen tvilsom.

4. Indeksverdi 0.62. Aggregeringsmetoden sier oss hvordan den GJENNOMSNITTLIGE endringen i breareal har vært for Norge sett under ett (jf metoden i NI).

Jeg tror alt. 2 er best, og mer riktig enn f.eks. alt. 4. Alt. 4 gir for eksempel Sørlandet en del vekt fordi det er en del fjell der. Men der er ingen breer der, hverken nå eller før. Aggregeringen bør bare skje over de arealene hvor indikatoren har dekning. På Sørlandet kan man si at indikatoren kunne hatt dekning om vi hadde noen grunn for å tro at brearealet økte, for da knune det ha økt inn i det eksisterende fjellarealet.

Indikatorverdien for Norge er da 0.7, som sier oss at 30% av brearealet er borte siden referanseperioden. Dette er med forbehold om at GØT er gitt fra en lineær skalering, men det kan godt tenkes at f.eks NVI vil ha grunnlag for å mene noe annet.

Usikkerhet

Jeg kan ikke se noen god måte å inkorporere usikkerhet inn i indikatorverdien. Det eneste jeg kan se for meg er en måte å lage en geografisk usikkerhet ved å først produsere indikatorverdier for mindre arealer (lage en raster med middels oppløning, f.eks. 10x10km2) og bootstrappe fra disse verdiene, men dette vil jeg ikke anbefale.